Code
library(tidyverse)
library(geomtextpath)
library(MoMAColors)
library(PNWColors)
library(bench)
library(patchwork)
library(stickylabeller)
sensible_sample <- function(x, ...){x[sample(length(x), ...)]}Simulation
Research question and summary of approach
Here, we are interested in whether changes in seasonal resource availability produce adaptive, plastic responses in the length of development of the sex role reversing Bush cricket K. nartee. Protandry, where males eclose as adults, or arrive to a breeding site earlier than females has been documented across numerous taxa, including several Bush cricket species (Simmons et al., 1994; Lehmann, 2012). However, it is unknown how the availability of resources during adolescence affects the degree of protandry observed in the population.
To generate testable predictions, we built an individual-based evolutionary simulation, where male and female eclosion time could evolve independently. We simulated breeding seasons in continuous time, but treated generations as non-overlapping annual events.
Within each generation, the simulation tracks a univoltine, gonochoristic population of adults that eclose into an environment where resources temporally increase in abundance. Females are immediately ready to breed upon eclosion, whereas males must first acquire resources, the energy from which they use to produce a spermatophore. As the breeding season progresses, males enter the mating pool and encounter receptive females at a certain rate, which leads to mating events. After mating, a female leaves the mating pool to oviposit fertilised eggs, the number of which depends on the resources currently available in the environment and the nutritional content of the spermatophore she consumed during mating. The female returns to the mating pool once all eggs from that reproductive bout have been laid. Mated males also leave the mating pool post mating, as they must produce a new spermatophore. Multiply-mated females store sperm from previous matings, producing sperm competition (see the Fitness section for details). Throughout the season, both sexes face a constant risk of mortality i.e. due to predation. Finally, the breeding season concludes abruptly, when all adults are killed by a heatwave, as is the case for this species in natural conditions.
Across generations, we then tracked the evolution of sex-specific emergence time. The composition of each generation is determined by individual reproductive success in the previous. We assume that the number of eggs laid in a season is far greater than the adult carrying capacity, which we constrain to be a constant \(N\), with all eggs equally likely to make it through density-dependent population regulation. Eggs therefore act as ‘lottery tickets’, with more fecund females and their mating partners holding more tickets (once again, see the Fitness section for details). We model emergence time as a quantitative trait determined by both an individual’s breeding value and environmental noise, with sex-limited expression. After \(G\) generations we recorded the mean trait value amongst females and males respectively, and quantified the level of protandry as the difference in mean breeding value between the sexes.
Importantly, we don’t model plasticity here. Instead, we allowed emergence to evolve towards the sex-specific optima in a given resource abundance scenario. Then, to test whether these optima change in response to resource availability, we ran a series of simulations, across which we varied the length of the breeding season. Specifically, we shifted the timing of the heatwave. At one extreme, the heatwave hits early in the season, when resources are still scarce. At the other, the heatwave hits very late, when ad libitum resources have been available for a considerable amount of time (akin to poor or high quality conditions during adolescence). If the environment experienced as a juvenile provides a cue for a sex-specific plastic response in emergence time, we expect this to shift sex-specific emergence times towards their phenotypic optima, as predicted by our simulations.
The remainder of this report details exactly how we built our individual-based model and includes all current results. It’s broken down into several sections (Emergence, Resource abundance, Fitness and running the simulation proper), within which I’ve tried to give a fuller description of the relevant process, while highlighting the biological assumptions we made to get the thing running.
Still to do
Thus far, I’ve only run simulations with a constant heatwave time across generations i.e. for a given simulation, the heatwave always hits at exactly time \(t\). A logical next step is to see what effect stochastic variation in the heatwave time has on the evolution of emergence time. Hanna and I suspect that bet-hedging will be important here.
I would also like to run a counterfactual simulation where males start in the mating pool, as I suspect this is why we never see protogyny.
Load packages
Emergence
This function turns a breeding value into a phenotypic value. We follow quant-gen convention and consider an individual’s phenotype as its breeding value, plus/minus some noise caused by the environment.
Code
emergence_sample <- function(emergence_breeding_value, p){
# T is the emergence breeding value
# we take advantage of the fact that the cumulative distribution of the emergence time is exp(t)/(exp(t)+exp(T)),
# hence the t that corresponds to a uniformly distributed p is ln(exp(T) p/(1-p)) = ln(exp(T))+ln(p/(1-p)) = T+ln(p/(1-p))
#p <- runif(length(emergence_breeding_value))
emergence_breeding_value+log(p/(1-p))
}Resource abundance
How resources increase with time is described by the logistic growth function
\[R(t) = \frac{1}{1 + e^{-k(t - t_0)}}\]
where \(t_0\) is the \(t\) value of the functions midpoint and \(k\) determines the rate at which resources enter the system. The function saturates at 1, which indicates the point where female fecundity is no longer limited by resource abundance. Here’s what it looks like for different values of \(t_0\), with \(k\) = 0.5.
Code
expand_grid(t = seq(from = -20, to = 20, by = 0.1),
t0 = seq(from = -15, to = 15, by = 2)) %>%
mutate(resources = 1/(1 + exp(-0.5*(t - t0))),
resources = case_when(t > 15 ~ NaN,
.default = resources)) %>%
ggplot(aes(x = t, y = resources)) +
geom_vline(xintercept = 15, linetype = 2, linewidth = 1) +
geom_line(aes(colour = t0, group = t0), linewidth = 1.5) +
scale_color_moma_c("Panton") +
scale_x_continuous(limits = c(-15, 16), expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = "Time, t",
y = "Resource abundance",
colour = "t0") +
theme_bw() +
theme(text = element_text(size = 16))
From this point on, I set \(k = 0.5\) and \(t_0 = 0\). I chose \(t_0 = 0\) because it’s also where I start mean emergence time for both sexes at generation zero. This just gives us a nice reference point for trait means as evolution progresses (neg trait means indicate evo of earlier emergence times). However time units are completely arbitrary and can easily be changed if this is not intuitive.
\(R(t)\) describes how resources flow into the system as time passes. Under the simplifying assumptions that spermatophore size is fixed across the population and that resource acquisition is independent of conspecific density and behaviour, the male refractory period is \(x-T\), where \(T\) is the time during the season when a male starts building a spermatophore and \(x\) is the time-point at which he has gathered sufficient resources to produce a complete spermatophore. We calculate spermatophore production between two points in the season as \[ A = \int_{T}^{x} \frac{1}{c_\mathrm{m} + e^{-k(t - t_0)}} \,dt \] where \(c_\mathrm{m}\) controls how efficiently males convert these resources into spermatophore. When \(c_\mathrm{m}\) is small, spermatophore production is efficient/cheap. We set the necessary amount of resources required to complete a spermatophore to \(A = 1\) and for a given value of \(T\) find males return to the mating pool at time
\[ x = \frac{kt_0 + \mathrm{Log}[\frac{-1+e^{c_\mathrm{m}k} + c_\mathrm{m}e^{c_\mathrm{m}k + k (T-t_0)}}{c_\mathrm{m}}]}{k} \]
Resource availability has the opposite effect on the female refractory period. We assume that females use the resources available to them immediately after mating to produce eggs, which they then spend time out of the mating pool ovipositing. Females lay eggs at a population-wide constant rate; the time spent out of the mating pool is therefore determined by the number of eggs they were able to produce immediately after mating. Longer refractory periods result from having a greater number of eggs to lay. Females also acquire resources from the spermatophore they receive during mating, which provides an additive ‘bump’ to the resources they have at their disposal. The efficiency at which females convert resources into eggs is controlled by the constant \(c_\mathrm{f}\), such that the female refractory is \(r_\mathrm{f}=c_\mathrm{f}(R(t) + s)\), where \(s\) is a population-wide constant that controls the nutritional content of the spermatophore. Any values of \(R(t) + s\) that exceed 1 are reset to 1; the value where female fecundity is no longer limited by resources.
An example
Set \(c_\mathrm{f}\) = 3, \(c_\mathrm{m}\) = 0.5, \(s\) = 0, \(k\) = 0.5, \(t_0\) = 0 and plot
Code
panel_1 <-
expand_grid(t = seq(from = -20, to = 20, by = 0.1),
t0 = 0,
k = 0.5) %>%
mutate(resources = 1/(1 + exp(-k*(t - t0))),
resources = case_when(t > 15 ~ NaN,
.default = resources),
end_y = 1/(1 + exp(-k*(-1 - t0)))) %>%
ggplot(aes(x = t, y = resources)) +
geom_segment(x = 0, xend=0,
y = -Inf, yend=0.5,
colour = "black", linetype = 2, linewidth = 0.7) +
geom_vline(xintercept = 15, linetype = 2, linewidth = 0.7) +
geom_line(linewidth = 1) +
scale_x_continuous(limits = c(-15.5, 15.5), expand = c(0, 0)) +
scale_y_continuous(expand = c(0.005, 0.005)) +
labs(x = "Time, t",
y = "Resource abundance") +
theme_bw() +
theme(text = element_text(size = 16),
legend.position = "none")
panel_2 <-
expand_grid(t = seq(from = -20, to = 20, by = 0.1),
t0 = seq(from = -15, to = 15, by = 1),
cm = 0.5,
fm = 3,
k = 0.5) %>%
mutate(male_timeout = ((k*t0 + log((-1 + exp(cm*k) + cm*exp(cm*k+k*(t-t0)))/(cm)))/k)-t,
female_timeout = (fm/(1 + exp(-k*(t - t0)))),
male_timeout = case_when(t > 15 ~ NaN,
.default = male_timeout),
female_timeout = case_when(t > 15 ~ NaN,
.default = female_timeout)) %>%
filter(t0 == 0) %>%
ggplot(aes(x = t)) +
geom_vline(xintercept = 15, linetype = 2, linewidth = 0.7) +
geom_vline(xintercept = 0, linetype = 2, linewidth = 0.7) +
geom_textline(aes(y = male_timeout, label = "Males"), size = 6, hjust = 0.12) +
geom_textline(aes(y = female_timeout, label = "Females"), size = 6, hjust = 0.12) +
scale_x_continuous(limits = c(-15.5, 15.5), expand = c(0, 0)) +
scale_y_continuous(expand = c(0.005, 0.005)) +
labs(x = "Time, t",
y = "Time out") +
scale_y_continuous(breaks = c(0, 2, 4, 6, 8, 10, 12,
14, 16, 18, 20)) +
theme_bw() +
theme(text = element_text(size = 16))
panel_1 / panel_2 + plot_layout(axes = "collect")
Build the refractory functions for the simulation
Code
# this doesn't find refractory period, but rather the time the male re-enters the mating pool
get_male_time_in <- function(c_m, t, t0, k){
(k*t0 + log((-1 + exp(c_m*k) + c_m*exp(c_m*k+k*(t-t0)))/(c_m)))/k
}
# females acquire resources via flowers and spermatophores
get_female_resource_quantity <- function(t, t0, k, s){
R_t <- (1/(1 + exp(-k*(t - t0)))) + s
if(R_t > 1){R_t <- 1}
R_t
}Fitness
For both sexes, the quantity maximised by selection is lifetime reproductive success. For females, this is equivalent to the total time they spend in time-out, ovipositing eggs. We treat eggs like tickets in a lottery, and draw mothers for adults in the next generation, weighted by the number of ‘tickets’ each female has. The relevant lottery is density-dependent viability selection, which we assume to occur prior to the reproductive life-stage, with each egg having equal probability of making it.
For males, reproductive success depends on mating success, the fecundity of their mates and the intensity of sperm competition. We model sperm competition under two different scenarios: 1) sperm-dumping, where a single mating is sufficient to fill the spermatheca and 2) complete sperm-mixing without any sperm dumping. In the sperm-dumping scenario, when a female mates with a second male, she dumps half of the total sperm she has received randomly (as she has received enough to fill her spermatheca twice-over), such that her spermatheca now contains equal parts of male one and male two’s ejaculate. If she mates a third time, she again dumps half the sperm, such that male three contributes half the present sperm, and male one and male two have their representation reduced to a quarter each. A female’s recent mates therefore have the highest chance of siring her offspring (there is last-male sperm precedence). To put this in perspective, a male’s probability of fathering offspring from a given reproductive bout reduces to ~0.01 if the female has subsequently mated six more times. 2) In the complete sperm-mixing scenario, a female’s spermatheca expands to store all the sperm she receives across her life. Paternity is evenly divided between her mates i.e. an egg of a thrice-mated female has a 33% chance of being fertilised by a sperm from each male. This is fair-raffle sperm competition. In both cases, I’ve made the simplifying assumption that fertilisation and sperm mortality have negligible effects on the number of sperm in the spermatheca. For completeness, we consider cases where females have mated up to 15 times, though this never occurs in the currently run simulations.
Code
# we model two different sperm competition scenarios
# 1) sperm dumping, leading to last male sperm precedence
# 2) fair raffle, where all inseminated sperm are stored by a female
# columns are how polyandrous a female has been
# rows hold the mating partner's paternity prob, which depends on how many times the female has subsequently mated
sperm_precedence_weights <-
data.frame(Mating1 = c(1, rep(0, 14)),
Mating2 = c(0.5, 0.5, rep(0, 13)),
Mating3 = c(0.5^2, 0.5^2, 0.5, rep(0, 12)),
Mating4 = c(0.5^3, 0.5^3, 0.5^2, 0.5, rep(0, 11)),
Mating5 = c(0.5^4, 0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 10)),
Mating6 = c(0.5^5, 0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 9)),
Mating7 = c(0.5^6, 0.5^6, 0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 8)),
Mating8 = c(0.5^7, 0.5^7, 0.5^6, 0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 7)),
Mating9 = c(0.5^8, 0.5^8, 0.5^7, 0.5^6, 0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5,
rep(0, 6)),
Mating10 = c(0.5^9, 0.5^9, 0.5^8, 0.5^7, 0.5^6, 0.5^5, 0.5^4, 0.5^3, 0.5^2,
0.5, rep(0, 5)),
Mating11 = c(0.5^10, 0.5^10, 0.5^9, 0.5^8, 0.5^7, 0.5^6, 0.5^5, 0.5^4, 0.5^3,
0.5^2, 0.5, rep(0, 4)),
Mating12 = c(0.5^11, 0.5^11, 0.5^10, 0.5^9, 0.5^8, 0.5^7, 0.5^6, 0.5^5,
0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 3)),
Mating13 = c(0.5^12, 0.5^12, 0.5^11, 0.5^10, 0.5^9, 0.5^8, 0.5^7, 0.5^6,
0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 2)),
Mating14 = c(0.5^13, 0.5^13, 0.5^12, 0.5^11, 0.5^10, 0.5^9, 0.5^8, 0.5^7,
0.5^6, 0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5, rep(0, 1)),
Mating15 = c(0.5^14, 0.5^14, 0.5^13, 0.5^12, 0.5^11, 0.5^10, 0.5^9, 0.5^8,
0.5^7, 0.5^6, 0.5^5, 0.5^4, 0.5^3, 0.5^2, 0.5))
fair_raffle_weights <-
data.frame(Mating1 = c(1, rep(0, 14)),
Mating2 = c(rep(1/2, 2), rep(0, 13)),
Mating3 = c(rep(1/3, 3), rep(0, 12)),
Mating4 = c(rep(1/4, 4), rep(0, 11)),
Mating5 = c(rep(1/5, 5), rep(0, 10)),
Mating6 = c(rep(1/6, 6), rep(0, 9)),
Mating7 = c(rep(1/7, 7), rep(0, 8)),
Mating8 = c(rep(1/8, 8), rep(0, 7)),
Mating9 = c(rep(1/9, 9), rep(0, 6)),
Mating10 = c(rep(1/10, 10), rep(0, 5)),
Mating11 = c(rep(1/11, 11), rep(0, 4)),
Mating12 = c(rep(1/12, 12), rep(0, 3)),
Mating13 = c(rep(1/13, 13), rep(0, 2)),
Mating14 = c(rep(1/14, 14), rep(0, 1)),
Mating15 = c(rep(1/15, 15)))
# mother function - mothers are drawn from the pool of egg-laying females, weighted by their relative fecundity
mother_finder_function <-
function(fecundity_matrix, # holds times for reproductive bouts
carrying_capacity){ # holds times when females end each egg-laying bout
# find mothers
# get likelihoods for mothering success
probs <- rowSums(fecundity_matrix, na.rm = T)/sum(fecundity_matrix, na.rm = T)
# find mothers of the next generation
mothers <- sensible_sample(x = 1:carrying_capacity, size = carrying_capacity, prob = probs, replace = T)
mothers
# sanity check - are all mothers female?
#if(setdiff(mothers, pop %>% filter(sex == 0) %>% pull(ID))>0){ # should return 0
# print("some mothers are male!")
#}
}
# father function
father_finder_function <-
function(mothers, # from mother finder function
fecundity_matrix,
mates,
sperm_comp_array){ # how does sperm comp work
fathers <- rep("NaN", length(mothers))
offspring <- 1 # pick a father for the first adult in the next gen; iterate through all individuals
# empty vector to be filled; females assumed to have mated a max of 15 times
fecundity_weights <- c(rep(0, nrow(sperm_comp_array)))
while(offspring <= length(mothers)){
# no. of mates
number_mates <- sum(!is.na(mates[mothers[offspring], ]))
# who are the mates
mate_IDs <- mates[mothers[offspring], 1:number_mates]
if(number_mates > 1){
fecundity_weights[1:number_mates] <-
fecundity_matrix[mothers[offspring], 1:number_mates] /
sum(fecundity_matrix[mothers[offspring], 1:number_mates])
paternity_weights <-
rowSums(data.frame(mapply(`*`,sperm_comp_array,fecundity_weights)))
father <- sensible_sample(mate_IDs, size = 1, prob = paternity_weights[1:number_mates])
}else{father <- mate_IDs}
fathers[offspring] <- father
offspring <- offspring + 1
}
fathers
}Define parameter space
Code
parameter_space <-
expand_grid(popsize = 1000, # 1000
mu_f_in = 0.1,
mu_f_out = 0.1,
mu_m_in = 0.1,
mu_m_out = 0.1,
v = c(8/(popsize / 2), 50/(popsize / 2)), # rate searching sex finds other sex
c_m = 0.5, # min time to regen spermatophore
c_f = seq(from = 0.5, to = 8, by = 0.5), # 0.5 to 8
s = c(0, 0.25, 0.5), # spermatophore boost 0, 0.25, 0.5
t0 = 0, # 50% resources at time 0
k = 0.5,
genome_wide_mutation = 0.2, # sd for normal with mu = 0 (variance is sigma^2)
heatwave_time = seq(from = 0, to = 15, by = 1), # when are adults killed each year
bet_hedging = c(0,1), # 0 = no, 1 = yes
gens = 2000)The complete simulation function
Code
simulation <- function(row, # which row of the parameter space
input, # the parameter space
sex_limited, # -1 = no, +1 = yes
OSR_recording_interval, # time points to record OSR
sperm_competition_weightings, # fair raffle, or last male precedence
male_start_out # TRUE or FALSE. Allows counter-factual test
){
#print(paste("doing row", row))
# define constants
popsize <- input$popsize[row]
mu_f_in <- input$mu_f_in[row] # female death rate in mating pool
mu_f_out <- input$mu_f_out[row] # female death rate out of the mating pool
mu_m_in <- input$mu_m_in[row] # male death rate in the mating pool
mu_m_out <- input$mu_m_out[row] # mating death rate out of the mating pool
v <- input$v[row] # velocity of male mate searching - all males have same value
c_m <- input$c_m[row]
c_f <- input$c_f[row]
s <- input$s[row]
t0 <- input$t0[row]
k <- input$k[row]
genome_wide_mutation <- input$genome_wide_mutation[row]
gens <- input$gens[row]
# constant heatwave time across generations
if(input$bet_hedging[row] == 0){
heatwave_time <- input$heatwave_time[row]
}
sperm_comp_weights <- sperm_competition_weightings
# initialise the population
# each generation fill the table with breeding values, emergence times and death times
if(sex_limited < 1){ # we only need to track one trait
population_attributes <-
data.frame(sex = rbinom(popsize, 1, prob = 0.5), # 0 is female, +1 is male
breeding_value = rnorm(popsize, mean = 0, sd = 1),
state = Inf,
p = runif(popsize)) %>%
mutate(emergence_time = emergence_sample(breeding_value, p)) %>% # get emergence time
arrange(emergence_time) %>% # order by emergence time to make downstream operations faster
mutate(ID = 1:popsize) # set ID for each individual
pop <-
population_attributes %>%
select(-c(p, breeding_value)) # remove un-needed columns
genetics <-
population_attributes %>%
select(ID, breeding_value)
} else{ # we need to track two traits
population_attributes <-
data.frame(sex = rbinom(popsize, 1, prob = 0.5), # 0 = female, +1 = male
f_limited_bv = rnorm(popsize, mean = 0, sd = 1),
m_limited_bv = rnorm(popsize, mean = 0, sd = 1),
state = Inf,
p = runif(popsize)) %>%
# get emergence time
mutate(emergence_time = case_when(sex > 0 ~ emergence_sample(m_limited_bv, p),
sex < 1 ~ emergence_sample(f_limited_bv, p))) %>%
arrange(emergence_time) %>% # order by emergence time to make downstream operations faster
mutate(ID = 1:popsize) # set ID for each individual
pop <-
population_attributes %>%
select(-c(contains("bv"), p)) # remove un-needed columns
genetics <-
population_attributes %>%
select(ID, f_limited_bv, m_limited_bv)
}
# for state:
# NaN: in mating pool,
# real number: time out with number indicating when they'll return
# Inf: dead or yet to emerge
# setup array to record results every generation
# first 4 cols track emergence means and sd,
# 5th is heatwave time
# 6th is prop females unmated in a generation
mean_trait_value <- matrix(NaN, ncol = 6, nrow = length(0:gens))
OSR_list <- list(rep(NaN, length(0:gens)))
# simulate evolution over many generations - within each generation, fitness is accrued during the breeding season, which we simulate with a Gillespie-like process of events
G <- 0
matings <- 1 # placeholder to get things started - sim will quit if no females mate
while(G <= gens & matings > 0){
print(paste("running generation", G))
# record mean breeding value in the population before selection
# genetic variance in trait values
if(sex_limited > 0){
mean_trait_value[G+1,1] <- mean(genetics$f_limited_bv) # f
mean_trait_value[G+1,2] <- sd(genetics$f_limited_bv) # f
mean_trait_value[G+1,3] <- mean(genetics$m_limited_bv) # m
mean_trait_value[G+1,4] <- sd(genetics$m_limited_bv) # m
} else{
mean_trait_value[G+1,1] <- mean(genetics$breeding_value) # both
mean_trait_value[G+1,2] <- sd(genetics$breeding_value) # both
}
# setup within gen necessities
#Mating_limit_exceeded <- "NO" # this will change to YES if a female mates more than 15 times, exceeding the limit we can track
# create an empty array to hold female mating times
# we assume that females can't mate more than 15 times across the season
female_out_times <- matrix(NaN, nrow = popsize, ncol = 15)
# and another that holds the male they mated with
female_mating_partners <- matrix(NaN, nrow = popsize, ncol = 15)
if(pop[1, 3] < -15){
t <- pop[1, 3] - 0.0001 # start tracking the population just before first emergence
} else{t <- -15} # start sim here at the latest
# set the heatwave time for the current gen
if(input$bet_hedging[row] > 0){
heatwave_time <- runif(1, min = input$heatwave_time[row], max =
15)
#max(input$heatwave_time)) #
mean_trait_value[G+1,5] <- heatwave_time
}
alive <- popsize # a stop early condition if everyone dies before the heatwave
# we record the OSR at regular intervals
next_OSR <- -15 # initial recording
OSR <- matrix(NaN, ncol = 2,
nrow = (OSR_recording_interval^-1)*length(next_OSR:heatwave_time))
OSR_row <- 1 # keep track of which row to update
#event_counter <- 0 # cut eventually
while(t < heatwave_time & # the timepoint when the heatwave kills all the adults
alive > 0){ # stop early if everyone dies before the heatwave
# gillespie through the events, tracking state changes
# find individuals out of the mating pool
# used to calculate next return to mating pool and death times
outfemales <- pop[pop$sex < 1 & is.finite(pop$state),]
outmales <- pop[pop$sex > 0 & is.finite(pop$state),]
# find individuals in the mating pool
# used to find the next potential mating and death times
infemales <- pop[pop$sex < 1 & is.na(pop$state),]
inmales <- pop[pop$sex > 0 & is.na(pop$state),]
# find next event
# emergence - remove past and select next one - pre-ordering makes this easy
next_emergence <- pop[pop$emergence_time > t, "emergence_time"][1]
# time in - Inf and NaN are possible options that the code can handle
next_time_in <- min(pop$state, na.rm = T)
# matings... if rate is 0, NaN produced.
pop_matingrate <- v*nrow(infemales)*nrow(inmales)
next_mating <- t + rexp(n = 1, rate = pop_matingrate)
# death... if rate is 0, NaN produced.
death_rate <-
mu_f_out*nrow(outfemales) +
mu_f_in*nrow(infemales) +
mu_m_out*nrow(outmales) +
mu_m_in*nrow(inmales)
next_death <- t + rexp(n = 1, rate = death_rate)
# find the next event and update t
t <- pmin(next_emergence,
next_time_in,
next_mating,
next_death,
next_OSR,
heatwave_time,
na.rm = TRUE) # ... if a rate is 0, NaN produced.
# record OSR
if(t == next_OSR & !is.na(next_OSR)){
OSR[OSR_row,1] <- nrow(infemales)
OSR[OSR_row,2] <- nrow(inmales)
next_OSR <- next_OSR + OSR_recording_interval # record OSR at specified time interval
OSR_row <- OSR_row + 1
}
# update the population by changing the states of individuals
# following emergence, initial state depends on sex
if(t == next_emergence & !is.na(next_emergence)){
# who is it
# because ID is ordered by emergence, select first row after filtering out old events
emergence_ID <- pop[pop$emergence_time >= t, "ID"][1]
if(pop[pop$ID == emergence_ID,]$sex > 0 & male_start_out){ # if male, find time male enters mating pool
pop[emergence_ID, "state"] <- get_male_time_in(c_m, t, t0, k)
} else{pop[emergence_ID, "state"] <- NaN} # if female or the counter-factual male case, start in time-in
}
# if the next event is a mating
if(t == next_mating & !is.na(next_mating)){
# see who mates
# everyone searches at the same rate
# so within sexes, all time-in individuals have the same probability of mating
# that means we can randomly sample one female and one male from the mating pool
who_mates_female <- sensible_sample(infemales$ID, 1)
who_mates_male <- sensible_sample(inmales$ID, 1)
# get_female_resource_quantity and calculate refractory period
R_f <- get_female_resource_quantity(t, t0, k, s) # use pre-defined function
f_refract_period <- R_f*c_f
# change state to time-out = real numbers
pop[who_mates_female, "state"] <- t + f_refract_period
pop[who_mates_male, "state"] <- get_male_time_in(c_m, t, t0, k)
# now update the female LRS table
# modify the recorded refractory period if egg-laying gets cut short by heatwave
if(f_refract_period + t > heatwave_time){
f_refract_period <- heatwave_time - t
}
# work out how many past repro bouts the newly mated female has had
# ID is the same as row number - females can thus be searched by row number
# for a specific row, we need to find the next column that has an NaN value
#e.g. if a female has mated 3 times previously, the next NaN will be in column 4
repro_bout <- which(is.nan(female_out_times[who_mates_female, ]))[1]
# Replace the next NaN in the row that corresponds to the appropriate female
female_out_times[who_mates_female, repro_bout] <- f_refract_period
female_mating_partners[who_mates_female, repro_bout] <- who_mates_male # record mate
}
# if the next event is a return to mating pool
if(t == next_time_in & !is.na(next_time_in)){
# find who it is and change state to time-in
pop[which.min(pop$state), "state"] <- NaN
}
# if the next event is a death, who dies
if(t == next_death & !is.na(next_death)){
who_dies <- NA # make sure this resets
out_f_weight <- (mu_f_out*nrow(outfemales))/death_rate
in_f_weight <- (mu_f_in*nrow(infemales))/death_rate
out_m_weight <- (mu_m_out*nrow(outmales))/death_rate
in_m_weight <- (mu_m_in*nrow(inmales))/death_rate
# which class does the death come from?
death_class <-
sensible_sample(c("out_f", "in_f", "out_m", "in_m"),
size = 1,
prob = c(out_f_weight, in_f_weight, out_m_weight, in_m_weight))
# all individuals within class have same death rate, choose 1 at random
if(death_class == "out_f"){
who_dies <- sensible_sample(outfemales$ID, 1)
# end the egg laying period prematurely
repro_bout_mu <- which(is.nan(female_out_times[who_dies, ]))[1] - 1 # find first NA, then go back a col
# subtract the time lost due to death from existing refractory period
if(pop$state[who_dies] > heatwave_time){
female_out_times[who_dies, repro_bout_mu] <-
female_out_times[who_dies, repro_bout_mu] - (heatwave_time - t)
} else{female_out_times[who_dies, repro_bout_mu] <-
female_out_times[who_dies, repro_bout_mu] - (pop$state[who_dies] - t)
}
}
if(death_class == "in_f"){
who_dies <- sensible_sample(infemales$ID, 1)}
if(death_class == "out_m"){
who_dies <- sensible_sample(outmales$ID, 1)}
if(death_class == "in_m"){
who_dies <- sensible_sample(inmales$ID, 1)}
# change state
pop[who_dies, "state"] <- Inf
}
# count individuals that haven't died
alive <- nrow(pop[pop$state != Inf | pop$emergence_time > t,])
#event_counter <- event_counter + 1
#print(paste("t =", round(t, 3), "event no. =", event_counter,
# "no. ind alive =", alive))
} # end within generation while loop
# once the season is complete, compute fecundity and male paternity
# from this we determine which breeding values are transmitted to the next generation
# length of time for each reproductive bout == fecundity
matings <- sum(!is.nan(female_out_times)) # if this is zero the sim will quit
if(matings > 0){ # if at least one female mated, produce the next generation
# find the mothers using our pre-written function
mothers <- mother_finder_function(female_out_times,
popsize)
# find the fathers using our pre-written function
fathers <- father_finder_function(mothers,
female_out_times,
female_mating_partners,
sperm_comp_weights)
# setup adults in the next generation
if(sex_limited > 0){ # two traits need to be tracked
next_gen <-
cbind(as.data.frame(mothers), as.data.frame(fathers)) %>%
merge(genetics, by.x = "mothers", by.y = "ID") %>%
merge(genetics, by.x = "fathers", by.y = "ID",
suffixes = c(".mother", ".father")) %>%
mutate(f_limited_bv = (f_limited_bv.mother + f_limited_bv.father)/2 +
rnorm(popsize, mean = 0, sd = genome_wide_mutation),
m_limited_bv = (m_limited_bv.mother + m_limited_bv.father)/2 +
rnorm(popsize, mean = 0, sd = genome_wide_mutation),
sex = rbinom(popsize, 1, prob = 0.5), # 0 is female, +1 is male
state = Inf,
p = runif(popsize),
# get emergence times
emergence_time =
case_when(sex > 0 ~ emergence_sample(m_limited_bv, p),
sex < 1 ~ emergence_sample(f_limited_bv, p))) %>%
arrange(emergence_time) %>% # order by emergence time to streamline later operations
mutate(ID = 1:popsize) # set ID for each individual
# get the two tables we need to run the next generation
pop <-
next_gen %>%
select(sex, state, emergence_time, ID) # remove un-needed columns
genetics <-
next_gen %>%
select(ID, f_limited_bv, m_limited_bv)
} else{ # one trait needs to be tracked
next_gen <-
cbind(as.data.frame(mothers), as.data.frame(fathers)) %>%
merge(genetics, by.x = "mothers", by.y = "ID") %>%
merge(genetics, by.x = "fathers", by.y = "ID",
suffixes = c(".mother", ".father")) %>%
mutate(breeding_value = (breeding_value.mother + breeding_value.father)/2 +
rnorm(popsize, mean = 0, sd = genome_wide_mutation),
sex = rbinom(popsize, 1, prob = 0.5), # 0 is female, +1 is male
state = Inf,
p = runif(popsize),
# get emergence times
emergence_time = emergence_sample(breeding_value, p)) %>%
arrange(emergence_time) %>% # order by emergence time to streamline later operations
mutate(ID = 1:popsize) # set ID for each individual
# get the two tables we need to run the next generation
pop <-
next_gen %>%
select(sex, state, emergence_time, ID) # remove un-needed columns
genetics <-
next_gen %>%
select(ID, breeding_value)
}
}
# record OSR across season
OSR_list[[G+1]] <- OSR
# record prop. unmated females
mean_trait_value[G+1,6] <- 1 - sum(!is.na(female_mating_partners[,1]))/(popsize - sum(pop[,1]))
# update the generation counter
G <- G + 1
# loop back to top of while loop...
}
results <- list(mean_trait_value, OSR_list)
results
}Run the simulation
Run the simulation for 10,000 generations, for a single set of parameters.
Code
parameters_single_example <-
expand_grid(popsize = 1000,
mu_f_in = 0.1,
mu_f_out = 0.1,
mu_m_in = 0.1,
mu_m_out = 0.1,
v = 8/(popsize / 2),
c_m = 1,
c_f = 8,
s = 0.25,
t0 = 0,
k = 0.5,
genome_wide_mutation = 0.2,
heatwave_time = 10,
bet_hedging = 0,
gens = 10000)
parameters_single_example_bh <-
expand_grid(popsize = 1000,
mu_f_in = 0.1,
mu_f_out = 0.1,
mu_m_in = 0.1,
mu_m_out = 0.1,
v = 8/(popsize / 2), # rate searching sex finds other sex
c_m = 1, # min time to regen spermatophore
c_f = 8,
s = 0.25, # spermatophore boost
t0 = 0, # 50% resources at time 0
k = 0.5,
genome_wide_mutation = 0.2, # sd for normal with mu = 0
heatwave_time = 10, # when are all the adults killed each year
bet_hedging = 1, # this is the only difference from above
gens = 10000)First let’s do it for our four cases, holding the heatwave time constant across generations
Code
# fair raffle
if(file.exists("results/single_run_fair_raffle.rds")){single_run_fair_raffle <- read_rds("results/single_run_fair_raffle.rds")
}else{single_run_fair_raffle <-
simulation(row = 1,
input = parameters_single_example,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = fair_raffle_weights,
male_start_out = TRUE)
write_rds(single_run_fair_raffle, "results/single_run_fair_raffle.rds")
}
# fair raffle with counter-factual
if(file.exists("results/single_run_fair_raffle_cf.rds")){single_run_fair_raffle_cf <- read_rds("results/single_run_fair_raffle_cf.rds")
}else{single_run_fair_raffle_cf <-
simulation(row = 1,
input = parameters_single_example,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = fair_raffle_weights,
male_start_out = FALSE)
write_rds(single_run_fair_raffle_cf, "results/single_run_fair_raffle_cf.rds")
}
# last male sperm precedence
if(file.exists("results/single_run_sperm_precedence.rds")){single_run_sperm_precedence <- read_rds("results/single_run_sperm_precedence.rds")
}else{single_run_sperm_precedence <-
simulation(row = 1,
input = parameters_single_example,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = sperm_precedence_weights,
male_start_out = TRUE)
write_rds(single_run_sperm_precedence, "results/single_run_sperm_precedence.rds")
}
# last male sperm precedence with counter-factual
if(file.exists("results/single_run_sperm_precedence_cf.rds")){single_run_sperm_precedence_cf <- read_rds("results/single_run_sperm_precedence_cf.rds")
}else{single_run_sperm_precedence_cf <-
simulation(row = 1,
input = parameters_single_example,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = sperm_precedence_weights,
male_start_out = FALSE)
write_rds(single_run_sperm_precedence_cf, "results/single_run_sperm_precedence_cf.rds")
}Now let’s do it for a mild bet-hedging case (uncertainty in heatwave time from time = 10 to time = 15), for our four scenarios
Code
# fair raffle
if(file.exists("results/single_run_fair_raffle_bh.rds")){single_run_fair_raffle_bh <-
read_rds("results/single_run_fair_raffle_bh.rds")
}else{single_run_fair_raffle_bh <-
simulation(row = 1,
input = parameters_single_example_bh,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = fair_raffle_weights,
male_start_out = TRUE)
write_rds(single_run_fair_raffle_bh, "results/single_run_fair_raffle_bh.rds")
}
# fair raffle with counter-factual
if(file.exists("results/single_run_fair_raffle_cf_bh.rds")){single_run_fair_raffle_cf_bh <-
read_rds("results/single_run_fair_raffle_cf_bh.rds")
}else{single_run_fair_raffle_cf_bh <-
simulation(row = 1,
input = parameters_single_example_bh,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = fair_raffle_weights,
male_start_out = FALSE)
write_rds(single_run_fair_raffle_cf_bh, "results/single_run_fair_raffle_cf_bh.rds")
}
# last male sperm precedence
if(file.exists("results/single_run_sperm_precedence_bh.rds")){single_run_sperm_precedence_bh <-
read_rds("results/single_run_sperm_precedence_bh.rds")
}else{single_run_sperm_precedence_bh <-
simulation(row = 1,
input = parameters_single_example_bh,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = sperm_precedence_weights,
male_start_out = TRUE)
write_rds(single_run_sperm_precedence_bh, "results/single_run_sperm_precedence_bh.rds")
}
# last male sperm precedence with counter-factual
if(file.exists("results/single_run_sperm_precedence_cf_bh.rds")){single_run_sperm_precedence_cf_bh <-
read_rds("results/single_run_sperm_precedence_cf_bh.rds")
}else{single_run_sperm_precedence_cf_bh <-
simulation(row = 1,
input = parameters_single_example_bh,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = sperm_precedence_weights,
male_start_out = FALSE)
write_rds(single_run_sperm_precedence_cf_bh, "results/single_run_sperm_precedence_cf_bh.rds")
}Run the simulation for 2,000 generations, for many combinations of parameters. To do this, we make use of the parallel package, which allows us to split the simulation effort across multiple cores.
For testing, I’ve first included code for running the simulation on a single core
Code
lapply(1:10, # rows to run
simulation,
input = parameter_space,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = sperm_mixing_weights)Run the proper simulation by multi-coring, splitting it into manageable chunks
Code
# split parameter space up
space_1 <- parameter_space[1:400,]
space_2 <- parameter_space[801:1600,]
#space_3 <- parameter_space[201:300,]
#space_4 <- parameter_space[301:400,]
#space_5 <- parameter_space[401:500,]
#space_6 <- parameter_space[501:600,]
#space_7_10 <- parameter_space[601:1000,]
#space_11_15 <- parameter_space[1001:1500,]
if(!file.exists("results/fair_raffle_1_400.rds")){
# get the number of cores on the computer
n.cores <- parallel::detectCores()
# make the cluster
cluster <- parallel::makeCluster(n.cores)
# split parameters across cluster and import required functions and libraries
parallel::clusterEvalQ(cluster, c(library(dplyr)))
parallel::clusterExport(cluster, c("space_1",
"emergence_sample",
"mother_finder_function",
"father_finder_function",
"get_female_resource_quantity",
"get_male_time_in",
"sensible_sample"))
results <- parallel::parLapply(cluster,
1:nrow(space_1),
simulation,
input = space_1,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = fair_raffle_weights)
write_rds(results, "results/fair_raffle_1_400.rds")
}
if(!file.exists("results/fair_raffle_801_1600.rds")){
# get the number of cores on the computer
n.cores <- parallel::detectCores()
# make the cluster
cluster <- parallel::makeCluster(n.cores)
# split parameters across cluster and import required functions and libraries
parallel::clusterEvalQ(cluster, c(library(dplyr)))
parallel::clusterExport(cluster, c("space_2",
"emergence_sample",
"mother_finder_function",
"father_finder_function",
"get_female_resource_quantity",
"get_male_time_in",
"sensible_sample"))
results <- parallel::parLapply(cluster,
1:nrow(space_2),
simulation,
input = space_2,
sex_limited = 1,
OSR_recording_interval = 0.5,
sperm_competition_weightings = fair_raffle_weights)
write_rds(results, "results/fair_raffle_801_1600.rds")
}Load the results
Code
results_1 <- read_rds("results/complete_results_1.rds")
results_2 <- read_rds("results/complete_results_2.rds")
results_3 <- read_rds("results/complete_results_3.rds")
results_4 <- read_rds("results/complete_results_4.rds")
results_5 <- read_rds("results/complete_results_5.rds")
results_6 <- read_rds("results/complete_results_6.rds")
results_7_10 <- read_rds("results/complete_results_7_10.rds")
results_11_15 <- read_rds("results/complete_results_11_15.rds")Plotting a single example
To get an idea for how long simulations need to run before trait values converge on some optimum, I ran a single simulation with the following parameter values:
Code
glimpse(parameters_single_example)Rows: 1
Columns: 15
$ popsize <dbl> 1000
$ mu_f_in <dbl> 0.1
$ mu_f_out <dbl> 0.1
$ mu_m_in <dbl> 0.1
$ mu_m_out <dbl> 0.1
$ v <dbl> 0.016
$ c_m <dbl> 1
$ c_f <dbl> 8
$ s <dbl> 0.25
$ t0 <dbl> 0
$ k <dbl> 0.5
$ genome_wide_mutation <dbl> 0.2
$ heatwave_time <dbl> 10
$ bet_hedging <dbl> 0
$ gens <dbl> 10000The evolution of protandry across many generations
Code
# wrangle data directly from simulations
plotting_data_fair_raffle <-
single_run_fair_raffle[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Fair raffle",
Male_start = "Begin in time out")
plotting_data_fair_raffle_cf <-
single_run_fair_raffle_cf[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Fair raffle",
Male_start = "Begin in time in")
plotting_data_sperm_precedence <-
single_run_sperm_precedence[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Sperm precedence",
Male_start = "Begin in time out")
plotting_data_sperm_precedence_cf <-
single_run_sperm_precedence_cf[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Sperm precedence",
Male_start = "Begin in time in")
# combine into single data frame
plotting_data_time_series <-
rbind(plotting_data_fair_raffle,
plotting_data_fair_raffle_cf,
plotting_data_sperm_precedence,
plotting_data_sperm_precedence_cf)
# make the plots
time_series_plots <-
plotting_data_time_series %>%
ggplot(aes(x = Generation,
y = mu,
group = Sex,
fill = Sex,
colour = Sex)) +
#geom_line() +
geom_ribbon(aes(ymin = mu - sd, ymax = mu + sd)) +
scale_fill_manual(values = c(pnw_palette("Sunset2", n = 5)[2],
pnw_palette("Sunset2", n = 5)[4])) +
scale_colour_manual(values = c(pnw_palette("Sunset2", n = 5)[2],
pnw_palette("Sunset2", n = 5)[4])) +
scale_x_continuous(limits = c(0,10000)) +
scale_y_continuous(limits = c(-7, 2),
breaks = c(-6, -4, -2, 0, 2)) +
labs(y = "Emergence time (standard deviation around mean breeding value)",
colour = "Sex") +
facet_wrap(Sperm_competition ~ Male_start, nrow = 2,
labeller = label_glue('Sperm comp type: {`Sperm_competition`}\nMales: {`Male_start`}')) +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
panel.border = element_rect(fill = NA, colour = "black", size = .8),
strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8))
matelessness_plots <-
plotting_data_time_series %>%
filter(Sex == "Female") %>%
ggplot(aes(x = Generation, y = Prop_females_unmated)) +
geom_line(alpha = 1) +
scale_x_continuous(limits = c(0,10000)) +
scale_y_continuous(limits = c(0,0.3)) +
labs(y = "Proportion females unmated") +
facet_wrap(Sperm_competition~Male_start) +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
panel.border = element_rect(fill = NA, colour = "black", size = .8),
strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8))
time_series_plots / matelessness_plots
Code
# wrangle data directly from simulations
plotting_data_fair_raffle_bh <-
single_run_fair_raffle_bh[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Fair raffle",
Male_start = "Begin in time out")
plotting_data_fair_raffle_cf_bh <-
single_run_fair_raffle_cf_bh[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Fair raffle",
Male_start = "Begin in time in")
plotting_data_sperm_precedence_bh <-
single_run_sperm_precedence_bh[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Sperm precedence",
Male_start = "Begin in time out")
plotting_data_sperm_precedence_cf_bh <-
single_run_sperm_precedence_cf_bh[1] %>%
as.data.frame() %>%
as_tibble() %>%
rename(Female_mu = X1,
Female_sd = X2,
Male_mu = X3,
Male_sd = X4,
Prop_females_unmated = X6) %>%
mutate(Generation = 0:10000) %>%
select(-X5) %>%
pivot_longer(cols = 1:4, names_to = "cat", values_to = "emergence_time") %>%
separate(cat, into = c("Sex", "Measure"), sep = "_") %>%
pivot_wider(names_from = Measure, values_from = emergence_time) %>%
mutate(Sperm_competition = "Sperm precedence",
Male_start = "Begin in time in")
# combine into single data frame
plotting_data_time_series_bh <-
rbind(plotting_data_fair_raffle_bh,
plotting_data_fair_raffle_cf_bh,
plotting_data_sperm_precedence_bh,
plotting_data_sperm_precedence_cf_bh)
# make the plots
time_series_plots_bh <-
plotting_data_time_series_bh %>%
ggplot(aes(x = Generation,
y = mu,
group = Sex,
fill = Sex,
colour = Sex)) +
#geom_line() +
geom_ribbon(aes(ymin = mu - sd, ymax = mu + sd)) +
scale_fill_manual(values = c(pnw_palette("Sunset2", n = 5)[2],
pnw_palette("Sunset2", n = 5)[4])) +
scale_colour_manual(values = c(pnw_palette("Sunset2", n = 5)[2],
pnw_palette("Sunset2", n = 5)[4])) +
scale_x_continuous(limits = c(0,10000)) +
scale_y_continuous(limits = c(-7, 2),
breaks = c(-6, -4, -2, 0, 2)) +
labs(y = "Emergence time (standard deviation around mean breeding value)",
colour = "Sex") +
facet_wrap(Sperm_competition ~ Male_start, nrow = 2,
labeller = label_glue('Sperm comp type: {`Sperm_competition`}\nMales: {`Male_start`}')) +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
panel.border = element_rect(fill = NA, colour = "black", size = .8),
strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8))
matelessness_plots_bh <-
plotting_data_time_series_bh %>%
filter(Sex == "Female") %>%
ggplot(aes(x = Generation, y = Prop_females_unmated)) +
geom_line(alpha = 1) +
scale_x_continuous(limits = c(0,10000)) +
scale_y_continuous(limits = c(0,0.3)) +
labs(y = "Proportion females unmated") +
facet_wrap(Sperm_competition~Male_start) +
theme_bw() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5),
panel.border = element_rect(fill = NA, colour = "black", size = .8),
strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8))
time_series_plots_bh / matelessness_plots_bh
The evolution of the OSR across many generations
For this example, we can also plot how the OSR changes in response to selection.
Code
OSR_data_fair_raffle <-
single_run_fair_raffle[2] %>%
as.data.frame() %>%
as_tibble() %>%
mutate(Time = seq(from = -15, to = 10.5, by = 0.5)) %>%
pivot_longer(cols = !last_col(),
names_to = "class",
values_to = "Individuals") %>%
separate(class, into = c("place", "generation"), sep = "\\.") %>%
mutate(sex = case_when(str_detect(place, "X1") ~ "Female",
.default = "Male"),
generation = if_else(is.na(generation), "0", generation),
generation = as.numeric(generation)) %>%
select(-place) %>%
filter(Time != 10.5) %>%
pivot_wider(names_from = sex, values_from = Individuals) %>%
mutate(OSR = Male / (Male + Female)) %>%
mutate(resources = 1/(1 + exp(-0.5*(Time - (0))))) %>%
filter(generation == 0 |
generation == 300 |
generation == 600 |
generation == 900 |
generation == 1200 |
generation == 1500 |
generation == 1800 |
generation == 2100 |
generation == 2400 |
generation == 3000 |
generation == 4000 |
generation == 5000 |
generation == 6000 |
generation == 7000 |
generation == 8000 |
generation == 9000 |
generation == 10000) %>%
mutate(Female = if_else(Female == 0, NaN, Female),
Male = if_else(Male == 0, NaN, Male),
Sperm_competition = "Fair raffle",
Male_start = "Begin in time out")
OSR_data_fair_raffle_cf <-
single_run_fair_raffle_cf[2] %>%
as.data.frame() %>%
as_tibble() %>%
mutate(Time = seq(from = -15, to = 10.5, by = 0.5)) %>%
pivot_longer(cols = !last_col(),
names_to = "class",
values_to = "Individuals") %>%
separate(class, into = c("place", "generation"), sep = "\\.") %>%
mutate(sex = case_when(str_detect(place, "X1") ~ "Female",
.default = "Male"),
generation = if_else(is.na(generation), "0", generation),
generation = as.numeric(generation)) %>%
select(-place) %>%
filter(Time != 10.5) %>%
pivot_wider(names_from = sex, values_from = Individuals) %>%
mutate(OSR = Male / (Male + Female)) %>%
mutate(resources = 1/(1 + exp(-0.5*(Time - (0))))) %>%
filter(generation == 0 |
generation == 300 |
generation == 600 |
generation == 900 |
generation == 1200 |
generation == 1500 |
generation == 1800 |
generation == 2100 |
generation == 2400 |
generation == 3000 |
generation == 4000 |
generation == 5000 |
generation == 6000 |
generation == 7000 |
generation == 8000 |
generation == 9000 |
generation == 10000) %>%
mutate(Female = if_else(Female == 0, NaN, Female),
Male = if_else(Male == 0, NaN, Male),
Sperm_competition = "Fair raffle",
Male_start = "Begin in time in")
OSR_data_sperm_precedence <-
single_run_sperm_precedence[2] %>%
as.data.frame() %>%
as_tibble() %>%
mutate(Time = seq(from = -15, to = 10.5, by = 0.5)) %>%
pivot_longer(cols = !last_col(),
names_to = "class",
values_to = "Individuals") %>%
separate(class, into = c("place", "generation"), sep = "\\.") %>%
mutate(sex = case_when(str_detect(place, "X1") ~ "Female",
.default = "Male"),
generation = if_else(is.na(generation), "0", generation),
generation = as.numeric(generation)) %>%
select(-place) %>%
filter(Time != 10.5) %>%
pivot_wider(names_from = sex, values_from = Individuals) %>%
mutate(OSR = Male / (Male + Female)) %>%
mutate(resources = 1/(1 + exp(-0.5*(Time - (0))))) %>%
filter(generation == 0 |
generation == 300 |
generation == 600 |
generation == 900 |
generation == 1200 |
generation == 1500 |
generation == 1800 |
generation == 2100 |
generation == 2400 |
generation == 3000 |
generation == 4000 |
generation == 5000 |
generation == 6000 |
generation == 7000 |
generation == 8000 |
generation == 9000 |
generation == 10000) %>%
mutate(Female = if_else(Female == 0, NaN, Female),
Male = if_else(Male == 0, NaN, Male),
Sperm_competition = "Sperm precedence",
Male_start = "Begin in time out")
OSR_data_sperm_precedence_cf <-
single_run_sperm_precedence_cf[2] %>%
as.data.frame() %>%
as_tibble() %>%
mutate(Time = seq(from = -15, to = 10.5, by = 0.5)) %>%
pivot_longer(cols = !last_col(),
names_to = "class",
values_to = "Individuals") %>%
separate(class, into = c("place", "generation"), sep = "\\.") %>%
mutate(sex = case_when(str_detect(place, "X1") ~ "Female",
.default = "Male"),
generation = if_else(is.na(generation), "0", generation),
generation = as.numeric(generation)) %>%
select(-place) %>%
filter(Time != 10.5) %>%
pivot_wider(names_from = sex, values_from = Individuals) %>%
mutate(OSR = Male / (Male + Female)) %>%
mutate(resources = 1/(1 + exp(-0.5*(Time - (0))))) %>%
filter(generation == 0 |
generation == 300 |
generation == 600 |
generation == 900 |
generation == 1200 |
generation == 1500 |
generation == 1800 |
generation == 2100 |
generation == 2400 |
generation == 3000 |
generation == 4000 |
generation == 5000 |
generation == 6000 |
generation == 7000 |
generation == 8000 |
generation == 9000 |
generation == 10000) %>%
mutate(Female = if_else(Female == 0, NaN, Female),
Male = if_else(Male == 0, NaN, Male),
Sperm_competition = "Sperm precedence",
Male_start = "Begin in time in")Code
pal_1 <- moma.colors("Alkalay2" , n=18, type="continuous")
OSR_f_fair_raffle <-
OSR_data_fair_raffle %>%
ggplot(aes(x = Time, y = Female,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of females in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_m_fair_raffle <-
OSR_data_fair_raffle %>%
ggplot(aes(x = Time, y = Male,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of males in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_fair_raffle <-
OSR_data_fair_raffle %>%
ggplot(aes(x = Time, y = OSR,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 0.01, width = 0) +
geom_line(alpha = 0.6, linewidth = 0.8) +
geom_line(aes(y = resources), colour = "black", linetype =2) +
labs(y = "OSR (prop. males in mating pool)",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
fair_raffle_OSR_plot <-
OSR_f_fair_raffle + OSR_m_fair_raffle + OSR_fair_raffle + guide_area() +
plot_layout(guides = 'collect')
fair_raffle_OSR_plot
Code
pal_1 <- moma.colors("Alkalay2" , n=18, type="continuous")
OSR_f_fair_raffle_cf <-
OSR_data_fair_raffle_cf %>%
ggplot(aes(x = Time, y = Female,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of females in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_m_fair_raffle_cf <-
OSR_data_fair_raffle_cf %>%
ggplot(aes(x = Time, y = Male,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of males in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_fair_raffle_cf <-
OSR_data_fair_raffle_cf %>%
ggplot(aes(x = Time, y = OSR,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 0.01, width = 0) +
geom_line(alpha = 0.6, linewidth = 0.8) +
geom_line(aes(y = resources), colour = "black", linetype =2) +
labs(y = "OSR (prop. males in mating pool)",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
fair_raffle_cf_OSR_plot <-
OSR_f_fair_raffle_cf + OSR_m_fair_raffle_cf + OSR_fair_raffle_cf + guide_area() +
plot_layout(guides = 'collect')
fair_raffle_cf_OSR_plot
Code
OSR_f_sperm_precedence <-
OSR_data_sperm_precedence %>%
ggplot(aes(x = Time, y = Female,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of females in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_m_sperm_precedence <-
OSR_data_sperm_precedence %>%
ggplot(aes(x = Time, y = Male,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of males in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_sperm_precedence <-
OSR_data_sperm_precedence %>%
ggplot(aes(x = Time, y = OSR,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 0.01, width = 0) +
geom_line(alpha = 0.6, linewidth = 0.8) +
geom_line(aes(y = resources), colour = "black", linetype =2) +
labs(y = "OSR (prop. males in mating pool)",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
sperm_precedence_OSR_plot <-
OSR_f_sperm_precedence + OSR_m_sperm_precedence + OSR_sperm_precedence +
guide_area() +
plot_layout(guides = 'collect')
sperm_precedence_OSR_plot
Code
OSR_f_sperm_precedence_cf <-
OSR_data_sperm_precedence_cf %>%
ggplot(aes(x = Time, y = Female,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of females in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_m_sperm_precedence_cf <-
OSR_data_sperm_precedence_cf %>%
ggplot(aes(x = Time, y = Male,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 1, width = 0)+
geom_line(alpha = 0.6, linewidth = 0.8) +
#geom_line(aes(y = resources), colour = "black", linetype =2) +
scale_y_continuous(limits = c(0, 300)) +
labs(y = "Number of males in mating pool",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
OSR_sperm_precedence_cf <-
OSR_data_sperm_precedence_cf %>%
ggplot(aes(x = Time, y = OSR,
group = generation, colour= generation)) +
geom_point(alpha = 1, height = 0.01, width = 0) +
geom_line(alpha = 0.6, linewidth = 0.8) +
geom_line(aes(y = resources), colour = "black", linetype =2) +
labs(y = "OSR (prop. males in mating pool)",
colour = "Generation") +
scale_colour_gradientn(colours = pal_1) +
theme_bw()
sperm_precedence_OSR_plot_cf <-
OSR_f_sperm_precedence_cf + OSR_m_sperm_precedence_cf + OSR_sperm_precedence_cf +
guide_area() +
plot_layout(guides = 'collect')
sperm_precedence_OSR_plot_cf
Full results
Load and wrangle the large output
Code
results_processer <-
function(row, results){
tibble(female_mean = results[[row]][[1]][1000,1],
male_mean = results[[row]][[1]][1000,3],
parameter_space_ID = row)
}
wrangled_results_1 <-
map_dfr(1:nrow(space_1), results_processer, results_1) %>%
bind_cols(space_1)
wrangled_results_2 <-
map_dfr(1:nrow(space_2), results_processer, results_2) %>%
bind_cols(space_2)
wrangled_results_3 <-
map_dfr(1:nrow(space_3), results_processer, results_3) %>%
bind_cols(space_3)
wrangled_results_4 <-
map_dfr(1:nrow(space_4), results_processer, results_4) %>%
bind_cols(space_4)
wrangled_results_5 <-
map_dfr(1:nrow(space_5), results_processer, results_5) %>%
bind_cols(space_5)
wrangled_results_6 <-
map_dfr(1:nrow(space_6), results_processer, results_6) %>%
bind_cols(space_6)
wrangled_results_7_10 <-
map_dfr(1:nrow(space_7_10), results_processer, results_7_10) %>%
bind_cols(space_7_10)
wrangled_results_11_15 <-
map_dfr(1:nrow(space_11_15), results_processer, results_11_15) %>%
bind_cols(space_11_15)
wrangled_results <-
bind_rows(wrangled_results_1,
wrangled_results_2,
wrangled_results_3,
wrangled_results_4,
wrangled_results_5,
wrangled_results_6,
wrangled_results_7_10,
wrangled_results_11_15)Build panels A-D
Code
# re-define the palette
pal_2 <- pnw_palette("Shuksan2",20)
heat_map_emergence <-
wrangled_results %>%
rename(female = female_mean, male = male_mean) %>%
pivot_longer(cols = female:male, names_to = "sex", values_to = "mean_emergence") %>%
filter(s != 0.1) %>%
ggplot(aes(x = c_f/c_m, y = heatwave_time)) +
geom_raster(aes(fill = mean_emergence)) +
scale_fill_gradientn(colours = pal_2) +
facet_wrap(sex ~ s, nrow = 2,
labeller = label_glue('Sex: {`sex`}\nSpermatophore value: {`s`}')) +
labs(x = "Female time-out / Male time-out\nwhen resources aren't limiting",
y = "Length of breeding season\n(heatwave time)",
fill = "Emergence time") +
scale_x_continuous(expand = c(0, 0),
breaks = c(3, 6, 9, 12, 15, 18)) +
scale_y_continuous(expand = c(0, 0)) +
theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
strip.text = element_text(size = 12),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.text = element_text(size = 12))Build panels E-F
Code
#pal_3 <- moma.colors("Alkalay2" , n=20, type="continuous")
pal_3 <- pnw_palette("Anemone",20)
heat_map_protandry <-
wrangled_results %>%
filter(s != 0.1) %>%
ggplot(aes(x = c_f/c_m, y = heatwave_time)) +
geom_raster(aes(fill = female_mean - male_mean)) +
scale_fill_gradientn(colours = rev(pal_3), limits=c(0,6)) +
facet_wrap(~ s, nrow = 1,
labeller = label_glue('Spermatophore value: {`s`}')) +
labs(x = "Female time-out / Male time-out\nwhen resources aren't limiting",
y = "Length of breeding season\n(heatwave time)",
fill = "Protandry") +
scale_x_continuous(expand = c(0, 0),
breaks = c(3, 6, 9, 12, 15, 18)) +
scale_y_continuous(expand = c(0, 0)) +
theme(panel.border = element_rect(fill = NA, colour = "black", size = .8),
strip.background = element_rect(colour = "black", fill = "Aliceblue", linewidth = .8),
strip.text = element_text(size = 12),
axis.title.x = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.text = element_text(size = 12))Plot
Code
heat_map_emergence/
heat_map_protandry + plot_layout(axes = "collect",
heights = c(2, 1))The left column of panels shows the results when the spermatophore has no nutritional value for females. This condition is clearly unrealistic (especially given males don’t re-enter the mating pool until they have produced a spermatophore) but it does provide an interesting comparison with the right column, where each spermatophore provides 25% of the resources females need to produce the maximum amount of eggs. At least for this comparison, the nutritional value of the spermatophore appears to have little qualitative impact on protandry. It’s possible that more nutritional spermatophores produce qualitative differences; that’s something I can explore.
For now though, let’s focus on the effects of two other parameters: the difference in time-out between the sexes when resources aren’t limiting (i.e. when \(R = 1\); x-axis) and the length of time resources are available before the heatwave hits (y-axis). For females (top row of panels), the results are intuitive. Emergence time depends on when the heat-wave hits: earlier heat-waves/late arriving resources = earlier relative emergence. In other words, if resources start appearing late, selection does not favour waiting for them to increase. For males, the effects are more complex (second row). Like females, males evolve earlier emergence times when the heatwave hits early, but the effect is more extreme. Additionally, if the difference in time-out with plentiful resources is small, male emergence is very sensitive to the timing of the heatwave. In this case, the extent of protandry (third row) is large with early heatwaves, but rapidly becomes minor when heatwaves occur later. If the difference in time-out is instead large, male emergence is much less sensitive to heatwave timing. The result is consistently strong protandry across a wide-range of heatwave timings.
Importantly, I don’t know if the ranges of the axis on these plots is realistic for the K. nartee system. What I’m calling small differences might in fact be large or vice versa.